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We perform several black-hole binary evolutions using fully nonlinear numerical relativity tech- 
niques at separations large enough that low-order post-Newtonian expansions are expected to be 
accurate. As a case study, we evolve an equal-mass nonspinning black-hole binary from a quasicir- 
cular orbit at an initial coordinate separation of D — lOOAf for three different resolutions. We find 
that the orbital period of this binary (in the numerical coordinates) is T — 6422M. The orbital 
motion agrees with post-Newtonian predictions to within 1%. Interestingly, we find that the time 
derivative of the coordinate separation is dominated by a purely gauge effect leading to an apparent 
contraction and expansion of the orbit at twice the orbital frequency. Based on these results, we 
improved our evolution techniques and studied a set of black hole binaries in quasi-circular orbits 
starting at D = 20M, D = 50M, and D = 100M for ~ 5, 3, and 2 orbits, respectively. We find 
good agreement between the numerical results and post-Newtonian predictions for the orbital fre- 
quency and radial decay rate, radiated energy and angular momentum, and waveform amplitude 
and phases. The results are relevant for the future computation of long-term waveforms to assist in 
the detection and analysis of gravitational waves by the next generation of detectors as well as the 
long-term simulations of black-hole binaries required to accurately model astrophysically realistic 
circumbinary accretion disks. 

PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION AND MOTIVATIONS 

Ever since the breakthroughs in Numerical Relativity 
(NR) of 2005 [IH3 , fully- nonlinear numerical simulations 
of merging black-hole binary (BHB) spacetimes have ad- 
vanced at a remarkable pace. From a simulation perspec- 
tive, BHB systems are characterized by the mass ratio 
(q = TO1/TO2), the individual spins of each BH (Si, 62), 
the initial separation (D), and the momenta of each BH 
(Pi, P2) (or equivalently the orbital eccentricity at some 
given separation). While all BHBs are computationally 
demanding, some of the interesting corners of the above 
parameter space remain particularly challenging. Initial 
investigations explored the orbital dynamics of similar- 
mass BHBs with low (or zero) spin, the computation of 
the associated gravitational radiation, and the character- 
ization of final remnant products of their mergers. These 
techniques were soon expanded to study mixed systems of 
black holes (BHs) and neutron stars 4-7], the merger of 
neutron-star binaries (forming BHs) [8- 10 , gravitational 
collapse [TT], and multiple black hole evolutions [T2] . 

Notable progress has recently been made in simulating 
binaries with mass ratios as small as q = 1/100 in [T3"lll4j . 
highly spinning black holes (with intrinsic spins S/m 2 > 
0.97) [T5], and ultra-relativistic collisions of black holes 
(reaching a Lorentz factor of 7 = 2.9) [16] . 

The close limit of black holes in numerical relativity 
has been studied in detail by the 'Lazarus' approach [171 — 
[23] . In this paper we perform a first study of BHBs 
at relatively "large" (from the relativistic point of view) 
initial separations. While we expect this problem to be 



well described by the post-Newtonian approximations to 
general relativity, it is important to assess how well the 
current numerical relativity techniques model this case, 
and to what extent they have to be adapted or corrected 
to reach a given accuracy. 

One of the main applications of the computation of 
accurate gravitational waveforms from BHBs is to as- 
sist gravitational wave observatories in detecting and an- 
alyzing gravitational radiation from astrophysical sys- 
tems; merging BHBs being a primary target [251 l2"5] - 
With the enhanced low-frequency sensitivity of Advanced 
LIGO compared to initial LIGO, much longer waveform 
templates are required as binaries will stay in the low- 
frequency part of the advanced LIGO sensitivity band for 
a much longer period. In addition, the more stringent ac- 
curacy criteria needed to match full numerical waveforms 
to post-Newtonian (or other analytical waveforms) [261 - 
[3~T] require up to a factor ten more waveform cycles than 
is the current norm. To explore the future possibility of 
generating such waveforms (and even longer) by purely 
full numerical techniques, we consider the case of a BHB 
at an initial separation of D = lOOAf . 

An independent motivation to study large-separation 
binaries comes from magnetohydrodynamics (MHD) 
studies of circumbinary disks around merging BHBs, 
which require many hundreds of orbital cycles before the 
accreting matter settles into a quasi-stationary state |32[ . 
BHB evolution can either provide a background space- 
time to evolve the MHD system, or the MHD and NR 
fields can be integrated in a self-consistent evolution. 

In Sec. |Hj we describe our numerical techniques to 
evolve black hole binaries using the moving puncture ap- 



proach [2J. In Sec. Ill we describe the results of evolving 
a BHB starting at a separation of D — 100M using the 
standard techniques developed for much closer binaries 
(D ~ 10M). We study the convergence of the orbital 
tracks of the holes for three resolutions (refined by suc- 
cessive factors of 1.2). We show results of new BHB sim- 
ulations of D = 100M , D = 50M, and D = 20M BHBs 
using improved evolution techniques in Sec. |IV| We also 
include in our analysis the proper separation between the 
two BHs, its decay, and the gravitational waveforms ex- 
traction at (at least) one wavelength from the sources. 
We compare all results with 3.5 post-Newtonian predic- 
tions. We conclude in Sec. [V] that current full numeri- 
cal techniques are able to accurately simulate the large 
separation regime but that the much longer time scales 
involved in the problem means that increases in simula- 
tions speeds by one to two orders of magnitude will be 
needed to complete the evolutions to merger. 



II. NUMERICAL SIMULATIONS 

We perform two different types of simulations here, a 
family of lower-accuracy, but high-speed, simulations of 
a BHB with initial separation of D = 100M at different 
resolutions, and a family of higher-accuracy simulations 
at a fixed resolution but with different initial separations. 
Note, here and below we use M to denote the unit of 
distance, mass, energy, etc., while we use m to denote the 
mass of an object. To avoid confusion between the mass 
of a BHB and the mass of its constituent BHs, we denote 
the masses of each BH in a binary with a subscript. 

We use the TwoPunctures thorn [33] to generate ini- 
tial puncture data [34] for the BHB simulations described 
below. These data are characterized by mass parame- 
ters m p (which are not the horizon masses), as well as 
the momentum and spin, of each BH. We evolve these 
BHB data sets using the LazEv [35] implementation of 
the moving puncture approach 2, 3 with the conformal 
function W = y/x = exp(— 2(f)) suggested by Ref. [56] . 
For the runs presented here, we use centered, eighth- 
order finite differencing in space 12J and a fourth-order 
Runge Kutta time integrator. (Note that we do not up- 
wind the advection terms.) Our code uses the Cac- 
tus/EinsteinToolkit [37H35] infrastructure. We use 
the Carpet [40] mesh refinement driver to provide a 
"moving boxes" style of mesh refinement. In this ap- 
proach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then 
moves these fine grids about the computational domain 
by following the trajectories of the two BHs. 

To reduce the computational costs, we performed an 
initial set of simulations (denoted by generation below) 
that use several low-accuracy approximations. Among 
them are the techniques introduced in Ref. [41] where 
the number of buffer zones at AMR boundaries is re- 
duced by lowering the order of finite differencing by suc- 
cessive orders near the AMR boundaries (here we use 6 



buffer zones), the use of simple interpolations of spec- 
tral initial data rather than using the complete spectral 
expansion 33J , and the copying the initial data to the 
two past time levels for use in prolongation at the initial 
timestep. All of these approximation proved to be useful 
for reducing the cost of numerical simulations, but each 
one also has the side-effect of introducing a (hopefully) 
small 0(h) error. 

For our higher-accuracy simulations (denoted by gen- 
eration 1 below), we use a full complement of 16 buffer 
zones at refinement boundaries, use the full TwoPunc- 
tures spectral expansion to generate initial data, and 
use the CARPETinit_3_timelevels option to populate 
the t < timelevels required for prolongation. We also 
used a strict 2:1 refinement in time, which is required for 
the init_3_timelevels option. In addition, we placed 
the outer boundaries at ~ 3000M, far enough to obtain 
waveforms. 

To compare the run speeds with and without the lower- 
order approximations, we performed two small runs using 
the grid structure and CFL factors of the generation 1 
simulations, one with the full complement of 16 buffer 
zones, and the other with only 6 buffer zones. The lat- 
ter was faster by 28% (136 M /day compared to 106 
M /day) The tests were performed on 20 nodes consist- 
ing of dual 6-core 3.47 GHz Nehalem CPUs. Increasing 
the CFL factor increases the run speed correspondingly 
(hence these low-accuracy techniques are quite useful in 
situations where high-accuracy is not required). 

We use a modified 1+log lapse and a modified T-driver 
shift condition [21I351I33], and an initial lapse a(t = 0) = 
2/(1 + if)% L ), where i^bl is the Brill-Lindquist conformal 
factor and is given by 
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where r, is the coordinate location of puncture i. The 
lapse and shift are evolved with 

{d t -p%)a = -2aK, (la) 

d t p a = (3/4)f a - a(i a , (lb) 

where we use a — 2 for the generation simulations 
presented below, while we use a spatially varying er = 
(o"o— Coo) exp(-r 4 /u; 4 )+cr 00 (where a = 2.0, a^ = 0.25, 
w = 80.0) for the generation 1 simulations. The reason 
for the differences is that, since we use strict 2:1 refine- 
ment in time for generation 1, we need a small a at large 
distances in order for the system to remain stable [44] . 

We use AHFinderDirect [JS] to locate apparent 
horizons. Here the spins of the two BHs are negligible 
and the horizon mass is equal to the irreducible mass 
(i.e. TO irr = yJA/{l&-K), where A is the surface area of 
the horizon). 



III. ZEROTH GENERATION RUNS 

Here we use our standard full numerical techniques (de- 
veloped to evolve BHB at initial separations D ~ 10M) 
to evolve a BHB at a separation of D = 100M and com- 
pare with the post-Newtonian predictions for the orbital 
trajectory. 



A. Quasicircular Post-Newtonian Analysis 

BHB quasicircular orbits at separations r ~ 100M 
(note that we use r to denote the post-Newtonian (PN) 
expression for the orbital separation) should be described 
accurately by the 2PN [46] expressions for the orbital fre- 
quency 
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and radial decay of the orbit (due to energy loss) 
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= mi + m-2 is the total mass and r\ = 
(toiTO2)/to 2 . The orbital period, T, at r = 100M for 
an equal-mass BHB is given by 
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6369M, 



(4) 



where mQ w 0.001, and m = \M. 

We thus estimate that the radial decay after one orbit 

is 



Ar = rT« -0.02M. 



(5) 



Note that in the above expressions r is the harmonic 
radial coordinate, which is related to the ADM-TT coor- 
dinate R (initially related to the numerical coordinates) 
by [5S] r = R + m 2 (2 + 29r))/(8R) + ■■■ . We also neglect 
the effects of the small eccentricities that are present in 
the full numerical runs. 



B. D — lOOAf Runs versus resolution 

In order to generate initial data, we need the punc- 
ture (coordinate) separation and momenta. We obtain 
approximations for the separation and momenta by us- 
ing the associated 3.5PN parameters. Our procedure is 
as follows: We start with 3.0 PN quasi-circular param- 
eters for an orbital separation of D ~ lOOAf . We then 
perform a 3.5PN evolution [47-45] of the orbital motion. 
The resulting orbit will be (moderately) eccentric. We 
then use the procedure in [50] to obtain successively bet- 
ter approximations of the initial momenta for this 3.5PN 
evolution until the eccentricity is reduced to the roundoff 



TABLE I: Initial data parameters for the full numerical sim- 
ulations. The punctures are located at ±(a;,0, 0), with mo- 
menta ±(p x ,p y ,0), spin (0,0,0), and puncture mass m p . run 
is the measured horizon masses, while T, is the period of the 
first orbit for the three resolution runs (0 for coarsest, 1 for 
medium, 2 for finest). mpN and Tpn are the corresponding 
post-Newtonian values 



a; = 50 

p x = -7.97939 x 10" 
p y = 2.55526 x 10~ 2 
m p = 0.49920645 
m H = 0.500615 



Tpn = 6365 
T = 6406 
Ti = 6420 
T 2 = 6422 
m PN = 0.500616 



TABLE II: The grid structure for the three numerical simu- 
lations presented here. The coarsest run had ft — fto = AM, 
while the two higher resolution runs had ft = fti — 4M/1.2 
and h = hi — 4M/1.44, respectively. 
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0.019 
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220 


ft/2 


0.038 
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110 
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0.076 
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55 
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0.152 


4 


25 


ft/16 


0.152 
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10 


ft/32 


0.304 
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ft/64 


0.304 


7 


2 


ft/128 


0.304 


8 


0.65 


ft/256 


0.304 



limit. After renormalizing to an ADM mass of 1M, we 
use the resulting PN parameters to solve for the initial 
data. The parameters are given in Table [I] 

We evolved the binary at three different resolutions 
characterized by a coarse-level resolution of fto = &M, 
fti = AM/1.2, and h 2 = M/f .44, respectively. For these 
runs, we choose a Courant-Friedrichs-Lewy (CFL) factor 
(i.e. the ratio dt/h between the timestep and gridsize) 
that varies with refinement level, as described in Table [H] 
In all cases, the boundaries were located at 400M and 9 
levels of (2:1) refinement were used. This means that the 
boundary is causally connected to the binary for most of 
the simulation. 

As seen in Fig. [T] the trajectory of the full numerical 
simulation is almost a closed circle with radius r ~ 100M. 
A more quantitative comparison is shown in Fig[2] where 
we plot the orbital frequency as a function of orbital sep- 
aration and compare to the 3.5PN predictions. As seen 
in Fig[2l the orbital frequency starts out much lower than 
predicted (which is purely a gauge effect) and then ap- 
proaches the PN value from below. 

A closer look at the time dependence of the coordinate 
separation (and its time derivative) shows the limited ac- 
curacy of the runs. A pure 3.5PN evolution indicates that 
this binary will take approximately t ~ 8.2 x 10 6 M to in- 
spiral to an orbital separation of r = 5M. During the first 
orbit, the 3.5PN binary inspirals only 0.02M from the ini- 
tial separation of r = 100M. As shown below, the full 
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FIG. 1: The orbital trajectory for the generation D — 100M 
simulations at different resolutions. 
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FIG. 2: The orbital frequency d<j> OI hit/dt versus orbital sep- 
aration r, as well as the Newtonian and 3.5 PN predictions 
for the generation D — 100M simulations. The mapping 
between r and Q is not unique because r is not a monotonic 
function of t. 



numerical binaries inspiral by two orders of magnitude 
more. While the 3.5PN simulation indicates that this bi- 
nary is not eccentric, there are significant oscillations in 
the orbital radius (usually associated with eccentricity). 
Interestingly, these oscillations occur at twice the orbital 
frequency (see Fig. |3| , while for eccentricity effects we ex- 
pect oscillations at the orbital frequency. In Figs. [3] and [IJ 
we plot the orbital separation and frequency versus time 
and compare to the PN prediction. Again, a significant 
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FIG. 3: The orbital separation versus time at three resolu- 
tions, as well as the 3.5PN prediction, for the generation 
D = 100M simulations. 
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FIG. 4: The orbital frequency versus time at three resolutions, 
as well as the 3.5PN prediction, for the generation D = 
100M simulations.. 



difference between the PN and numerical results is appar- 
ent. We will show in the next section that this behavior 
is mostly a coordinate artifact and proper distance mea- 
surements greatly ameliorate this issue. The coordinate 
artifact may be produced by the zero speed gauge mode 
(131 contained in our choice of the shift condition ( lb I . 



In a previous study |51j , we found that unphysical in- 
creases or decreases in the horizon mass can have a signif- 
icant impact on the orbital trajectory. In Fig. [5] we plot 
the horizon mass as a function of resolution versus time. 
We observe an increase in mass, which may be the cause 
of the artificial inspiral. The mass increase is small, but 
this secular growth (which gets worse as the resolution 
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FIG. 5: The horizon mass of the individual BHs versus time 
of the generation simulations for the three resolutions. Both 
the secular growth with time for all resolution, and the growth 
in mass as the gridspacing is reduced, indicate a loss of accu- 
racy at later times for these generation simulations. 



is increased) can affect the accuracy of longer term runs. 
We will see how to better control this artifact in the next 
section. 

Finally, as seen in Fig. [61 the convergence order, as 
measured by comparing the differences in the orbital sep- 
arations and phase between the ho and hi runs and the 
hi and hi runs, appears to be higher than eighth order. 
In the past, when we encountered this apparent super 
convergence, we found in the end that the error was an 
oscillatory function of the grid resolution (see [55]). How- 
ever, if we attempt to extrapolate, assuming eighth-order 
convergence, we find that the error in the orbital separa- 
tion increases to about 0.34Af and then appears to level 
off. The phase error increases to about 0.025 rad. 



IV. FIRST GENERATION RUNS 

In our study of waveform accuracy for closer BHBs 52J , 
we found that the low-accuracy approximations used in 
our generation simulations may have a larger impact 
than expected. As seen in Fig. [5j mass gain is an issue 
for generation 0. By correctly initializing the grid and by 
lowering the CFL factor to 0.25, we mitigate this mass 
loss issue. Importantly, we also place the outer bound- 
aries a factor of ~ 10 farther away. Note that it is the 
CFL factor in the inner regions (near each BH) that is 
responsible for the conservation (or lack thereof) of the 
BH mass. So even though the generation 1 runs have 
a larger CFL factor in the outer zones then the genera- 
tion runs, they are still more accurate in terms of mass 
preservation. 




FIG. 6: (Left) Convergence plots for the orbital phase 
(top-left) and separation (bottom-left) showing stronger than 
eighth-order convergence. If a quantity F(h) is eighth-order 
convergent, then F(h ) - F(fti) = 4.29982[F(fti) - F(h 2 )]. 
Here we see hyperconvergence, which almost certainly means 
that the errors is oscillating with h rather than behaving 
monotonically. The top-right plot shows the phase error of the 
coarsest resolution run (assuming eighth-order convergence 
and comparing the Richardson extrapolated value), while the 
bottom-right plot shows the orbital separation error. 



As an additional analysis tool, we use the proper dis- 
tance between the two horizons as measured along the 
coordinate line joining the two punctures |53) . which we 
call the semiproper distance, or SPD, below (note that 
the minimal geodesic does not follow this line). 

For the generation 1 runs, we performed simulations 
with a grid structure similar to the low resolution ho 
generation run, but with several additional coarser lev- 
els, and studied BHBs with initial orbital separations of 
D = 20M, D = 50M, and D = 100M. We placed the 
outer boundaries at 3200M and used 14 levels of refine- 
ment with a coarsest resolution of h = 32M (extending 
the grid to 6400Af stably can be achieved by halving the 
gauge baseline CToo or extending the grid without refining 
it by a factor 2 or halving the CFL condition). We used 
strict 2:1 refinement in time and a full complement of 16 
buffer zones at the boundary of each refinement level. 

The conservation of the BH masses during the evo- 
lution is an important indicator of the accuracy of the 
full numerical runs, and is particularly important for 
long term evolutions. We show the horizon mass for the 
D = lOOAf generation 1 run and compare with the gener- 
ation runs in Fig. [7] We observe that while generation 
runs show an strong growth of the masses with time 
and with increasing resolution, the profile of the genera- 
tion 1 runs is rather flat, and we expect this to be still 
the case for higher resolutions as observed in Ref. |52) . 
Fig. 2. 



A. Varying the BHB initial separation 

On of the robust features of both our generation and 
generation 1 runs is the approach of the full numerical 
orbital frequency to the PN predicted values from be- 
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FIG. 7: Comparing the D = 100M mass conservation for the 
generation and generation 1 runs. Note that the generation 
1 simulation shows no late time growth. The conservation of 
the mass at levels of below 1 part in 10 5 is crucial for long 
term evolutions. 



TABLE III: The eccentricity and orbital period of the numer- 
ical simulations, the 3 PN quasicircular orbital periods at the 
initial separation and separation after one orbit, and the 3.5 
PN orbital inspiral after one orbit. 



D/M 



Tnr/M 



Tpn/M 



Ar/M 



20 


0.0003 


615 


50 


0.0002 


2313 


100 


0.0006 


6422 



601-592 

2283-2279 

6370-6368 



-0.208 
-0.054 
-0.02 



low (after the numerical coordinate settle into a quasi- 
stationary regime [which takes ~ 1 orbit]). For the 
D = 20M case (see Fig. [8]), the numerical and PN predic- 
tions agree to two significant digits (e.g., SIpn = 0.01142, 
while fi num = 0.01125, at D = 18.84M) after the gauge 
settles down. The progression is similar for the cases of 
orbits starting at separations D = 50M and D — 100M, 
where the simulations completed 2.5 and 1.5 orbits, re- 
spectively. The retrograde features in these figures are 
due to the remaining (small) eccentricity of the initial 
quasicircular data and gauge effects. We report in Ta- 
ble |lll] the eccentricities of these BHBs, as measured via 
the formula e s ~ s 2 s, where s(t) is the semiproper dis- 
tance, as well as the measured orbital period and the 
3 PN prediction, including the decay rate Ar after one 
orbital period. 

As we saw in the generation runs (see Fig. pi), the 
coordinate separation of the holes shows an interesting 
gauge effect, namely an oscillation in time at twice the 
expected frequency. In order to show that this is a non- 
physical effect, we computed the proper distance along 
the line joining the black holes, i.e., the SPD measure 
described above (we use this measure since an accurate 
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FIG. 8: The orbital frequency d(j> OI bit / dt versus orbital sep- 
aration r and versus shifted semiproper distance s, as well 
as the Newtonian and 3.5 PN predictions for the D = 
20M, 50M, 100M generation 1 runs. Here s~— s - (s - r ), 
where so is the initial SPD and ro is the initial PN orbital 
separation. 



minimal proper distance between horizons is computa- 
tionally very intensive). The results are shown in Fig. [9] 
We first observe that the proper distance is much closer 
to being monotonic and featureless than the coordinate 
separation. We also see that there are residual oscil- 
lations (of much lower amplitude than those observed 
in the coordinate separation) at roughly the orbital fre- 
quency, which likely are due to the eccentricity, as well 
as a possible small residual gauge effect. 

In Fig. [9| we subtract a constant offset from the proper 
distances in order to make comparison with the coordi- 
nate separation easier. We see the striking agreement 
(apart from the overall constant offset between the two) 
at later times, i.e. after the first orbit, showing that the 
coordinate gauge effects settle down after one cycle, quite 
independent from the initial radius of the orbit, i.e. in 
all cases studied, D/M = 20, 50, fOO. 

In Fig. |10| we compare the orbital frequency, as a 
function of time, for the generation and generation 1 
(D = lOOAf ) simulations. Note how the generation 1 ap- 
pears to outperform the highest resolution generation 
run. 

Using our more accurate generation 1 simulations, and 
the SPD measure of the orbital separation, we study in 
more detail the orbital decay of the numerical runs and 
compare them with the PN predictions. Here were ex- 
amine s(t) as a function of time, where s(i) is the SPD. 
In Fig. [Tl] we see the good agreement between the mean 
values of s(t) and the 3.5PN predictions. The agreement 
is slightly better for the closer D/M = 20 case, which we 
evolved for three orbits, than for the more demanding 
D/M = 50, 100 cases, which did not complete two or- 
bits. There is an oscillation corresponding closely to the 
orbital period of the D/M = 20, 50, 100 runs (according 
to PN, T/M = 598, 2281, 6369, respectively). Note that 
in Fig. [Tl] we used the same ordinate scale and there- 
fore the amplitudes of oscillations are comparable, which 
indicates that the eccentricities of the farther separated 
BHBs are larger. 

In order to have a visual notion of how tight the orbits 
look in proper distance s space, we plot the semiproper 
trajectory r{4>) — s(cos</>, sin</>). Here <fi is the coordi- 
nate azimuthal angle, while s is the SPD. Our results 
are shown in Fig. [l2j The scale of the zoom-in plots is 
the same for all runs (making comparison of the rela- 
tive tightness of the orbits straightforward). Note that 
the larger separation between successive orbits for the 
D = 100M case than for the D = 50M case likely reflects 
some combination of residual gauge effects and trunca- 
tion errors (see Fig. 111. While the D — 100M orbital 



D=20M 



decay rate appears to be too large by a factor of 2-3, the 
gauge invariant radiated energy is consistent with PN 
predictions (see Sec. IV B and Fig. 111. 



We note that, according to 3.5PN evolutions, the D = 
20M BHB completes 36 orbits and takes t = 14, 200M to 
inspiral to D = 5M. The D = 50M BHB, on the other 
hand completes 370 orbits and takes t — 531,716M to 
get to the same separation. Finally, the D = 100M BHB 
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FIG. 9: Proper distance (minus a constant) and coordinate 
distance comparison for the D = 20M, 50Af , 100M genera- 
tion 1 runs. The D = lOOAf plot also shows the coordinate 
separation versus time for the generation runs. 
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FIG. 10: The orbital frequency for the generation and gener- 
ation 1 runs versus time (D — 100A/) and the PN prediction. 
Note how the generation 1 appears to outperform the highest 
resolution generation run. 



completes 2064 orbits and takes t = 8, 223, 170M . 

B. Comparison with Post-Newtonian waveforms 

The leading (2, 2)-mode of the PN waveform [M] is 
given by 

robs ft(a,2) = ^mJ^rjxH^e- 2 ^ (6) 

where ip = (p — 2iV/il ln(Sl/Slo) , f is the orbital phase, 

H^ = l+(- 1 °l + 55 r 1 ] x + 27rx 3 / 2 
\ 42 42 V 
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1512 

irri - 24ii] ) x 5/2 



(7) 



x = (mil) 2 ' 3 , and firj is an arbitrary positive constant 
(see [54]) which we took to be the initial orbital frequency. 
For quasicircular orbits, up to 3 PN, the orbital frequency 
is given by [57] 
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where R denotes the ADM-TT radial coordinate. 
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FIG. II: Simple proper distance decay rate versus post- 
Newtonian prediction for the D - 2QM,D = 50M, and 
D — 100M generation 1 runs. Note that the periods of os- 
cillations roughly match those of the orbital period, T/M = 
598, 2281, 6369 respectively 
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FIG. 12: The orbital trajectory obtained using the coordi- 
nate orbital phase and the proper orbital separation along the 
coordinate line separating the two BHs for the generation 1 
simulations. The zoomed plots show the tightness of the spi- 
ral. 



For large r, ^4 is given by [55J 

?/>4 = \ (h+ - ihy) « -2m 2 fl 2 h {2 > 2) 



(9) 



The radiated energy and angular momentum are give 
by Eqs. (7.4a)- (7.4b) of Ref. [55] 
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for quasicircular orbits. Hence, we 



can approximate the energy E and angular momentum 
J radiated per quasicircular orbit by 



dt 



dE dE2n 

dt dt Q' l J 



with a similar expression for the radiated angular mo- 
mentum (where the integral is over the time interval cor- 
responding to two complete cycles of the {I = 2,m = 2) 
mode of ipi, and hence one orbit of the BHB). 

Comparisons between full numerical and PN wave- 
forms can provide a strong test of both approaches. For 
the generation 1 runs, we extended the computational 
domain to a cube of sides 6400M in order to be able 
to extract the gravitational waveforms in the radiation 
zone. The effects of the extraction radii on the ampli- 
tude and phase of the numerical waveforms is apparent 
in Fig. 13 where we plot the D = 20M the waveforms as 
seen by observers at R ohs /M = 100, 300, 400, 800. Note 
that in this case, the half orbital period (measuring the 
period of the waves) is ~ 300Af . We see that extracting 
at R bs/M = 100 leads to inaccuracies in both the am- 
plitude and phase, but for R ohs /M = 300, 400, 800 the 
waveforms line up; indicating that extracting at least at 
one wavelength from the system is necessary. In all cases, 
we used the perturbative extrapolation formula given in 
Ref. [57]. 

lim[r^ m M)] 

r— >-oo 



-o(iCbs) 



r=flo bs 

(12) 



where for our choice to the tetrad we have ^4 = (1/2 — 
M/r) tpf um . For larger orbital separations the orbital pe- 
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FIG. 13: The waveform extracted at different observer loca- 
tions for the D = 20M generation 1 run. Here the waveforms 
have been extrapolated to oo using Eq. (12 I. While we see 
that the extraction at R b s = 100 has large phase errors, the 
phases and amplitudes stabilize for -Robs > 300, i.e. one wave- 
length distance from the sources. 




4000 



FIG. 14: Comparison of ripi and the Null variable ^ for the 
D — 50AI simulation. The boxes indicate regions where the 
fast Fourier transform smoothing operation on the CCE data 
distorted the waveform. 



riod increases (See Eqs. (2M)) approximately as ~ D 3 / 2 
hence we need observer locations of at least R b s /M = 
1150 and R b s /M = 3200, respectively, for binaries with 
separations D/M 
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as 



50 and D/M = 100. In Fig 
an independent validation of the finite radius extraction, 
we compare the perturbative extraction as defined above 



(12 1 with the Cauchy characteristic extraction (CCE) 
58] 



code described in Ref. [58] and observe the good agree- 
ment among them in the common region of validity. 

We find a striking agreement (particularly at larger 
separations) between these full numerical waveforms and 
the PN waveforms given by Eq. fcty. In Fig. 15 we plot 



TABLE IV: Energy and angular momentum radiated per or- 
bit (initial and second). 



D 



20M 
50M 
100M 



D 



20M 
50M 
100M 



Sm/Mrmm 

(5.68 ±0.02) x 1(T 
(2.4 ±0.1) x 10~ 6 
(2.3 ±0.4) x 10~ 7 



Sm/M-p-N 



5.43-5.62 x 10 
2.52-2.53 x 10 
2.36-2.36 x 10 



SJ/Ml 



OJ/iKI num 

(5.39 ±0.01) x 10" 3 



JTJWn 



V J.9±0.1) x 10" 

(2.4 ±0.4) x 10" 



5.20-5.30 x 10 
9.16-9.18 x 10 
2.39-2.39 x 10 



Rip4 extrapolated using Eq. (12) and the 3.5PN [55] pre- 



diction. The differences in the phase and amplitudes are 
within the numerical noise. For the D = 20 case, we 
still observe improvements between the 2PN expression, 
as truncated in Eq. (|7|, and the 3.5PN expression. At 
the larger initial separations, D = 50M, 100M, even the 
lower-order PN expressions [i.e. Eq. |7|] show excellent 
agreement with the full numerical waveform. 

We note that based on the generation convergence 
simulations, we may expect that the orbital phase error 
in the D = lOOAf simulation run is under 0.025 rad (the 
argument being that the generation 1 simulations should 
be more accurate at a given resolution than the gener- 
ation simulations). Based on the observed agreement 
between the PN and NR simulation, it appears that the 
phase error is indeed (relatively) small. 

Finally, we give the energy and angular momentum ra- 
diated, per orbit, in Table |TV) For the table, we measure 
the energy and angular momentum radiated over one pe- 
riod of the (£ = 2,rn = 2) mode (however, we use all 
modes up to I — 4 to calculate the energy and angular 
momentum radiated), we then multiply by 2 (we do this 
because we do not have two full cycles at i? h s = 3000Af 
for D = 100M). We approximate the error in these mea- 
surements by calculating 5m and 5 J at i? h s = 3000M 
(R ohs = 1600M for D = 20M) and R ohs = 1500M 
(R ohs = 800M for D = 20M) and take the difference 
between these two measurements as the error. We com- 
pute the energy and angular momentum radiated at the 
initial separation and after one orbit including the radial 
decay given by Ar w ( ^^ ) T, where 
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V. DISCUSSION 

We performed a first exploration of full numerical evo- 
lutions of a black-hole binaries with large initial sepa- 



11 



D=20M 



2e-05 



1e-05 - 



-1e-05 - 



-2e-05 



5e-07 



2.5e-07 



-2.5e-07 



-5e-07 




1000 2000 3000 4000 5000 

t/M 

D=100M 



4e-08 



2e-08 



- 



-2e-08 - 



-4e-08 




2000 4000 6000 8000 10000 

t/M 



FIG. 15: Comparison of 3.5PN and full numerical waveforms 
for the D = 20M, 50M, 100M generation 1 simulations. Note 
the excellent agreement in both amplitude and phase and the 
very small scale of the amplitudes. 



rations in order to evaluate how numerical techniques 
developed for close binaries (initial separations of around 
R ~ 10M) handle this regime. We studied a prototypical 
binary with an initial separation of D = 100M. Given 
the large orbital period involved (T s=s 6400M), we re- 
stricted our evolutions to the first few orbits. We find 
that the full numerical simulations agree with the post- 
Newtonian predictions for the gravitational waveform, 
orbital frequencies, orbital decay rate, and radiated en- 
ergy. These are nontrivial features given the length scales 
involved in the problem and the very small amplitude of 
the radiation. 

From these first studies we can draw several conclu- 
sions: 

i) The initial pulse of spurious radiation (and gauge re- 
laxation) still requires a period of the order of one orbital 
cycle to settle (that seems to be quite independent of the 
initial orbital radii). This implies longer evolution times 
to obtain accurate waveform information. Alternatively, 
one could use initial data with some information of the 
realistic radiation content along the lines of Ref. [59 . 
ii) Given the long evolution times required and the long 
wavelengths involved, the location of the computational 
boundaries should allow for the extraction of radiation 
at (at least) one wavelength from the sources. We note 
that more efficient techniques to treat the evolution in 
this far zone include the use of more accurate boundary 
conditions [50], multi-patch schemes [HI], and the choice 
of coordinates that are better adapted to the problem 

m 

iii) A pure 3.5PN evolution indicates that this binary 
will take approximately t ~ 8.2 x 10 6 Af to inspiral to 
an orbital separation of D — 5AI. During this inspiral, 
the binary would complete 2064 orbits. Our full numeri- 
cal simulations on 20 nodes (3.46GHz dual Intel proces- 
sors with 6 cores each) produced an average evolution 
of nearly 100M per day. This would led to a total of 
over 200 years to complete the evolution. Dramatic im- 
provements in the speed of the evolution codes, possibly 
using use hardware accelerators [62] , or novel numerical 
techniques, such as implicit-explicit methods [53], will be 
needed in order to make simulations from these separa- 
tions to mergers possible. 

iv) Based on the generation 1 results and our results 
in |52j . we can foresee a generation 2 set of runs that 
use the current runs to reduce eccentricity using for in- 
stance the method |50] , have the computational bound- 
aries moved to even larger radii (which would require 
setting (Too to a smaller value), use higher-order AMR 
prolongation and numerical dissipation, and replace the 
semiproper distance with the proper distance of the short- 
est geodesic joining the two horizons. 

We note that the eccentricity reduction method of [50] , 
when applied to s(t) for the D = 100M configuration, 
gives a very small change in the initial tangential mo- 
mentum and a very large change in the initial radial mo- 
mentum. In particular, we find Sp t /pt — 6.056 x 10~ 6 
and Sp r /p r = 18.6. The latter result is surprising as 
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it is not consistent with the magnitude of the PN radial 
momentum. We speculate that this indicates that the os- 
cillations in s(t) are not due to purely eccentricity effects, 
but have a significant gauge component, as well. If one 
were to modify the PN momentum as predicted here, and 
use these data as the starting point for a 3.5PN evolution, 
the resulting binary would have an initial eccentricity of 
5.5 x 10 , which is roughly the eccentricity we measure 
for the numerical binary. 

In conclusion we have shown that we can use cur- 
rent numerical techniques to accurately model the quasi- 
adiabatic evolution of black-hole binaries at radii of the 
order of D = 100M and generate the gravitational wave- 
forms from these binaries, but the speed of the time inte- 
gration techniques will need to be improved by two orders 
of magnitude before such simulations can be routine. 
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